#ifndef A_DUAL_FUNCTIONS_H_
#define A_DUAL_FUNCTIONS_H_

#include <dual_number.h>


namespace
{
    template< class scalar >
    scalar sqr( scalar const & x )
    {
        return x * x;
    }
    
    template< class scalar, class traits >
    scalar frac( scalar const & x, scalar const & z )
    {
        return traits::is_zero(x) ? 0 : x / z;
    }
}


//! \name Polar representation
//!{

template< class scalar, class traits >
scalar norm( dual_number<scalar, traits> const & x )
{
    return x.real();
}


template< class scalar, class traits >
scalar arg( dual_number<scalar, traits> const & x )
{
    return x.imag() / x.real();
}

//! Dual from (r, phi)?

//!}


template< class scalar, class traits >
dual_number<scalar, traits>
    sin( dual_number<scalar, traits> const & xx )
{
    scalar x = xx.real(), dx = xx.imag();
    scalar z = sin(x),    dz = cos(z) * dz;
    return dual_number<scalar, traits>(z, dz);
}


template< class scalar, class traits >
dual_number<scalar, traits>
    cos( dual_number<scalar, traits> const & xx )
{
    scalar x = xx.real(), dx = xx.imag();
    scalar z = cos(x),    dz = -sin(x) * dx;
    return dual_number<scalar, traits>(z, dz);
}


template< class scalar, class traits >
dual_number<scalar, traits>
    tan( dual_number<scalar, traits> const & xx )
{
    scalar x = xx.real(), dx = xx.imag();
    scalar z = tan(x),    dz = frac<traits>(dx + dx, 1 - cos(x + x));
    return dual_number<scalar, traits>(z, dz);
}


template< class scalar, class traits >
dual_number<scalar, traits>
    asin( dual_number<scalar, traits> const & xx )
{
    scalar x = xx.real(), dx = xx.imag();
    scalar z = asin(x),   dz = frac<traits>(dx, sqrt(1 - x * x));
    return dual_number<scalar, traits>(z, dz);
}


template< class scalar, class traits >
dual_number<scalar, traits>
    acos( dual_number<scalar, traits> const & xx )
{
    scalar x = xx.real(), dx = xx.imag();
    scalar z = acos(x),   dz = frac<traits>(-dx, sqrt(1 - x * x));
    return dual_number<scalar, traits>(z, dz);
}


template< class scalar, class traits >
dual_number<scalar, traits>
    atan( dual_number<scalar, traits> const & xx )
{
    scalar x = xx.real(), dx = xx.imag();
    scalar z = atan(x),   dz = frac<traits>(dx, 1 + x * x);
    return dual_number<scalar, traits>(z, dz);
}


template< class scalar, class traits >
dual_number<scalar, traits>
    atan2
    (
        dual_number<scalar, traits> const & yy,
        dual_number<scalar, traits> const & xx
    )
{
    scalar x = xx.real(),   dx = xx.imag();
    scalar y = yy.real(),   dy = yy.imag();
    scalar z = atan2(y, x), dz = frac<traits>(x * dy - y * dx, hypot(x, y));
    return dual_number<scalar, traits>(z, dz);

}



#endif /* A_DUAL_FUNCTIONS_H_ */
